package edu.gmu.atelier;

/**
 * Performs the DFT and IDFT in matrix form. 
 * I'd say that each signal frequency we know is a vector in the matrix.
 * The the x vector determines how much of each signal goes into the
 * output signal (i.e. the signal that is sampled), Wx=s.  
 * <p>
 *                [x0]   [s0]
 * [w0 w1 ... wn] [x1] = [s1]
 *                [x2]   [s2]
 *                [xn]   [sn]
 * <p>
 * Of course, we are given the sampled vector s=(s0,s1,...,sn)^T and 
 * need to compute the weighting coefficient (amplitude) 
 * vector x=(x0,x1,...,xn)^T. We do not need to decompose - the ws
 * are orthonormal and therefore W^-1 = W^T, x = (W^T)s.
 * <p>
 * Looked at at this way, after multiplying (W^T)s, the x vector
 * contains the amplitude determined for each corresponding signal
 * frequency (q1 == s1, q2 == s2,...).
 * <p>
 * Each row "signal" is "correlated" using inner product to produce
 * the amplitude at that signal frequency.
 * <p>
 *                  [s0]   [w0^T s] = [w01 s0 + w02 s2 + w0n sn] = [x0]
 * [w0 w1 ... wn]^T [s1] = [w1^T s] = [w11 s0 + w12 s2 + w1n sn] = [x1]
 *                  [s2]   [w2^T s] = [w21 s0 + w22 s2 + w2n sn] = [x2]
 *                  [sn]   [w3^T s] = [w31 s0 + w32 s2 + w3n sn] = [x3]
 * <p>
 * Really the only difficulty is creating the signal frequency matrix W.
 * <p>
 * Interestingly, I believe there is a connection with least-squares 
 * projections - is vector s within the column space of Wx?  Intuitively
 * I'd say the answer is NO!  So why do we always get a result?
 * <p>
 * Note W^T W x = W^T s is a projection of s onto the 
 * column space of W (?) Need to confirm, may not be "correct"
 * because we can exactly reproduce s from x implying that s was 
 * within the span(W) (projection would only be able to 
 * reproduce s_hat, closest point to s). Therefore, my conclusion
 * is that span(W) fills all of C^n in which vector s is in also, 
 * which was evident by the fact that W consists of n independent
 * basis vectors.
 * 
 * @author James H. Pope
 */
public class DFTSolver //implements Solver
{
    private Matrix w = null;
    private Matrix wt= null;
    
    private int fs = 0;
    private int n = 0;
    
    public DFTSolver( int n, int fs )
    {
        this.fs = fs;
        this.n  = n;
        //double intervalSpacing = fs / (double)n;
        /*
         * Generate n different signals evenly spaced over fs
         * The signals are orthogonal to each other, if we make
         * them orthonormal, can easily determine inverse W^-1=W^T
         * Note that since W consists of complex vectors, the transpose
         * also involves taking the conjugate, this is handled within
         * the complex vector implementation but realize it happens
         * so more precisely W^-1=conj(W^T)=W^H.
         */
        ComplexVector[] ws = new ComplexVector[n];
        for (int k = 0; k < ws.length; k++)
        {
            ComplexElement[] wi = new ComplexElement[n];
            // In some way frequency is captured within this formulation
            double basisSignal = (2.0 * k * Math.PI) / n;
            /*
             * Now generate representative sample
             * values for signal at different times
             */
            for (int t = 0; t < wi.length; t++)
            {
                // Moved (2.0 * k * Math.PI)/n to outer loop
                //double ck = (2.0 * t * k * Math.PI) / n;
                double ck = basisSignal * t;
                double real =    Math.cos(ck);
                double imag = -1*Math.sin(ck); // why -1? already conjugate?
                wi[t] = new ComplexElement(real, imag);
            }
            ws[k] = new ComplexVector(wi);
            //ws[k] = new ComplexVector(wi).normalize();
            // I believe the 1/n scales to make normal?
            System.out.println("Norm "+k+" is "+ws[k].norm2());
        }
        this.w = new Matrix(ws, ComplexVectorFactory.INSTANCE);
        this.wt= this.w.transpose();
    }
    
    public Matrix getW()
    {
        return this.w;
    }
    
    public Matrix getWT()
    {
        return this.wt;
    }
    
    //------------------------------------------------------------------------//
    // Interface methods
    //------------------------------------------------------------------------//
    public Vector solveFreq( Vector s ) // a.k.a. DFT
    {
        // That was easy :)
        // Yes, we need to normalize!
        //return b.mult(wt);
        return s.mult(wt).mult( new ComplexElement(1.0/n, 0) );
    }
    
    public Vector solveTime( Vector x ) // a.k.a. IDFT
    {
        // That was easy :)
        // No need to scale by n?
        //return x.mult(w).mult( new ComplexElement(n, 0) );
        return x.mult(w);
    }
    
}
